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Abstract.  Two  different  hybrid  particle-continuum  methods  are  described  for  simulation  of  nonequilibrium  gas  and 
plasma  dynamics.  The  first  technique,  used  for  nonequilibrium  hypersonic  gas  flows,  uses  either  a  continuum 
description  or  a  particle  method  throughout  a  flow  domain  based  on  local  conditions.  This  technique  is  successful  in 
reproducing  the  results  of  full  particle  simulations  at  a  small  fraction  of  the  cost.  The  second  method  uses  a 
continuum  model  of  the  electrons  combined  with  a  particle  description  of  the  ions  and  atoms  for  simulating  plasma 
jets.  The  physical  accuracy  of  the  method  is  assessed  through  comparisons  with  data  measured  in  space.  These 
examples  illustrate  that  the  multi-scale  physical  phenomena  associated  with  nonequilibrium  conditions  can  be 
simulated  with  physical  accuracy  and  numerical  efficiency  using  such  hybrid  approaches. 
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INTRODUCTION 

Nonequilibrium  conditions  occur  in  the  gas  and  plasma  flows  of  many  real  systems.  One  form  of 
nonequilibrium  concerns  the  generation  of  non-Maxwellian  velocity  distribution  functions  (VDFs).  Physically, 
such  nonequilibrium  is  produced  by  very  strong  gradients  in  flow  field  properties,  for  example  ion  shock  waves 
and  boundary  layers,  and  by  rarefied  flow  conditions.  Another  form  of  nonequilibrium  concerns  different 
species  in  the  gas  or  plasma  having  very  different  VDFs  from  one  another.  In  this  article,  a  summary  is  provided 
of  on-going  development  and  application  of  two  different  hybrid  fluid-particle  methods  for  simulating  these 
types  of  nonequilibrium  flows.  The  first  method  combines  solution  of  the  Navier-Stokes  equations  using 
traditional  methods  from  Computational  Fluid  Dynamics  (CFD)  with  the  particle-based  direct  simulation  Monte 
Carlo  (DSMC)  technique  [1].  This  hybrid  method  decides,  based  on  local  flow  conditions,  whether  to  use  CFD 
or  DSMC  in  each  region  of  a  flow  field.  The  DSMC  technique  is  more  physically  accurate  in  strongly 
nonequilibrium  regions,  but  is  at  least  an  order  of  magnitude  slower  than  CFD.  The  performance  of  this  method 
in  terms  of  physical  accuracy  and  numerical  efficiency  is  illustrated  through  its  application  to  hypersonic  gas 
flows.  The  second  hybrid  method  uses  fluid  and  particle  methods  everywhere  to  simulate  the  plasma  jets  created 
by  spacecraft  propulsion  systems  [2] .  In  this  technique,  electrons  are  modeled  using  a  fluid  approach  in  which 
their  conservation  equations  are  solved  using  CFD  methods.  The  heavy  species  (ions  and  atoms)  are  modeled 
using  a  combination  of  two  particle  methods.  The  DSMC  technique  is  again  employed  to  simulate  collisions 
while  the  Particle  In  Cell  (PIC)  method  is  used  to  self-consistently  accelerate  ions  in  electro- static  fields.  The 
performance  of  this  hybrid  method  is  illustrated  through  its  application  to  simulate  plasma  jets  created  by  a 
spacecraft  Hall  thruster. 


HYBRID  METHOD  FOR  NONEQUILIBRIUM  GAS  FLOWS 

For  relatively  low  altitude  portions  of  a  hypersonic  vehicle  entry  trajectory,  the  atmospheric  density  is 
relatively  high,  and  simulations  of  flows  around  hypersonic  vehicles  should  be  performed  using  traditional 
Computational  Fluid  Dynamics  (CFD)  by  solving  either  the  Euler  or  the  Navier-Stokes  (NS)  equations.  At  very 
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high  altitudes,  at  the  edge  of  an  atmosphere,  the  density  is  low  such  that  there  are  very  few  collisions  between  the 
molecules  and  atoms  in  the  flow  around  the  vehicle.  This  rarefied  flow  regime  can  be  computed  using  the  direct 
simulation  Monte  Carlo  (DSMC)  method  [3].  Relatively  speaking,  CFD  methods  for  solving  the  NS  equations 
are  about  an  order  of  magnitude  faster  than  the  DSMC  method.  However,  the  lack  of  collisions  in  the  low 
density  gas  makes  the  physics  of  the  NS  equations  invalid. 

A  flow  configuration  that  is  characterized  as  being  in  the  continuum  regime  overall  can  often  contain 
localized  regions  of  rarefied  flow.  Hypersonic  flow  over  a  cylinder  provides  an  excellent  example  of  this 
behavior,  and  the  relevant  phenomena  are  illustrated  in  Fig.  2.  In  such  a  flow,  localized  regions  of  rarefied  flow 
may  occur  in:  (1)  the  bow  shock  wave;  (2)  the  boundary  layer;  and  (3)  the  wake.  In  the  cases  of  the  bow  shock 
and  boundary  layer,  the  local  rarefaction  is  caused  by  the  very  strong  gradients  that  produce  small  localized 
length  scales  that  in  turn  produce  large  localized  Knudsen  numbers.  In  the  wake,  the  density  is  a  small  fraction 
of  that  in  the  forebody  flow  and  rarefaction  is  caused  here  directly  by  the  associated  large  mean  free  path. 

To  compute  such  flows  with  existing  computational  tools  we  face  the  dilemma  of  either:  (1)  using  CFD  and 
accepting  that  there  will  be  some  (unknown)  error  associated  with  inaccurate  physical  modeling  of  the  locally 
rarefied  regions;  or  (2)  trying  to  use  the  DSMC  technique  and  accepting  the  incredibly  high  computational  cost, 
especially  for  three  dimensional  flows.  Thus,  either  CFD  or  DSMC  on  its  own  fails  to  provide  a  comprehensive 
computational  modeling  capability  across  all  flow  regimes  encountered  by  a  hypersonic  vehicle.  A  natural 
solution  to  this  problem  is  to  develop  a  hybrid  simulation  technique  that  employs  a  CFD  method  for  as  much  of 
the  flow  field  as  possible  (due  to  its  superior  numerical  performance)  that  switches  to  using  DSMC  only  in  those 
regions  of  the  flow  where  the  physics  description  provided  by  the  CFD  method  is  inadequate. 

Continuum  Computational  Fluid  Dynamics 

Under  continuum  flow  conditions,  the  traditional  methods  of  Computational  Fluid  Dynamics  (CFD)  can  be 
employed  in  which  the  basic  conservation  equations  are  solved  numerically.  LeMANS  is  a  hypersonic  CFD  code 
that  solves  the  2D/3D  Navier-Stokes  equations  using  a  line  implicit  method  on  general  unstructured  meshes.  The 
code  is  parallelized  using  domain  decomposition.  Thermo-chemical  non-equilibrium  effects  are  included  by 
solving  a  separate  finite-rate  vibrational  energy  equation,  as  well  as  individual  species  conservation  equations 
that  include  source  terms  due  to  finite-rate  chemistry.  The  ability  to  simulate  a  weakly  ionized  gas  is  included. 
Further  details  of  the  code  can  be  found  in  [4] . 

The  Particle-Based  Direct  Simulation  Monte  Carlo  Method 

The  DSMC  technique  uses  the  motions  and  collisions  of  particles  to  perform  a  direct  simulation  of 
nonequilibrium  gas  dynamics  [3].  A  key  step  in  DSMC  is  the  separation  of  particle  motion  and  collision.  The 
particles  are  first  moved  through  physical  space  by  the  product  of  their  individual  velocities  and  a  time-step  that 
is  smaller  than  the  local  mean  free  time.  After  movement,  the  particle  locations  are  fixed,  the  particles  are 
collected  into  cells  that  have  dimensions  of  the  local  mean  free  path,  and  in  each  of  these  cells  a  number  of 
collisions  is  processed  consistent  with  local  flow  conditions.  Macroscopic  flow  properties  such  as  density  and 
temperature  are  obtained  by  time- averaging  particle  properties  in  the  cells  over  several  thousand  iterations.  The 
DSMC  technique  is  the  standard  numerical  method  for  simulating  rarefied  hypersonic  flows. 

MONACO  is  a  parallel  implementation  of  the  DSMC  method  in  which  a  computational  cell  is  taken  as  the 
basic  unit  of  the  simulation  [5].  Using  the  cell  as  the  basic  unit,  rather  than  using  the  particle  as  in  the  usual 
DSMC  algorithm,  provides  great  flexibility  in  terms  of  using  unstructured  grids  and  parallel  domain 
decomposition.  MONACO  contains  models  for  a  variety  of  physical  phenomena  including  collisional 
momentum  exchange,  rotational  energy  exchange,  vibrational  energy  exchange,  chemical  reactions,  and  wall 
collisions. 


A  Hybrid  CFD-DSMC  Method 

A  hybrid  CFD-DSMC  approach  developed  at  the  University  of  Michigan  is  called  the  Modular  Particle 
Continuum  (MPC)  method  [1]  and  merges  existing  CFD  and  DSMC  codes.  The  accuracy  of  a  hybrid  particle- 
continuum  method  relies  on  the  proper  positioning  of  the  DSMC-CFD  interfaces.  The  interface  must  lie  in  near¬ 
equilibrium  regions  where  solution  of  the  NS  equations  will  introduce  minimal  or  no  error.  Typically,  particle 


and  continuum  regions  are  determined  by  evaluating  a  continuum  breakdown  parameter  in  the  flow  field.  The 
MPC  method  uses  the  gradient-length  Knudsen  number: 

^nGLL~~Q^Q  (1) 

where  Q  represents  local  flow  quantities  such  as  density,  temperature,  or  velocity  magnitude,  X  is  the  local  mean- 
free-path,  and  the  gradient  is  evaluated  in  all  directions  to  determine  the  maximum  gradient.  It  has  been  shown 
for  hypersonic  flows  [6]  that,  in  regions  of  the  flow  field  where  KnGLL  <0.05,  the  discrepancy  between  a  NS  and 
DSMC  solution  is  less  than  5%.  Thus,  these  regions  could  be  solved  using  a  NS  solver  with  little  error.  The  MPC 
method  begins  with  a  NS  solution  of  the  entire  flow  field  and  then  uses  Eq.  (1)  to  decompose  the  domain  into 
CFD  and  DSMC  regions. 

The  MPC  method  uses  state-based  coupling  to  transfer  information  between  particle  and  continuum  regions. 
After  evaluation  of  the  breakdown  parameter,  the  particle  region  is  extended  by  a  few  extra  cells  into  the 
continuum  domain  to  create  an  overlap  region.  The  overlap  regions  help  to  improve  accuracy  of  the  DSMC 
fluxes.  Next,  one  row  of  NS  and  two  rows  of  DSMC  boundary  cells  are  initialized.  The  DSMC  and  NS  domains 
are  then  coupled  by  transferring  information  across  the  interfaces.  With  state-based  coupling,  this  involves 
updating  the  boundary  conditions  of  each  solver.  In  this  way,  information  transfer  into  both  the  particle  and 
continuum  regions  is  handled  through  existing  boundary  procedures  already  used  by  both  solvers.  The  DSMC 
boundary  cells  are  continually  filled  with  particles  consistent  with  the  flow  properties  in  the  corresponding  NS 
cell  using  the  Chapman-Enskog  distribution  based  on  the  local  macroscopic  state  and  gradients,  known  from  the 
NS  solver.  As  particles  in  the  DSMC  domain  interact  and  their  distributions  evolve  in  time,  the  MPC  method 
also  tracks  the  macroscopic  variation  in  each  DSMC  cell.  In  order  to  provide  these  averaged  properties  with 
minimal  statistical  scatter,  a  mixture  of  spatial  and  temporal  averaging  is  used.  Specifically,  MPC  uses  the  sub¬ 
relaxation  technique  proposed  by  [7].  These  averaged  DSMC  properties  are  then  used  to  update  boundary 
conditions  for  the  NS  solver. 

Figures  1  illustrate  the  performance  of  the  hybrid  code  for  computation  of  normal  shock  waves  of  argon. 
Normalized  density  profiles  for  a  Mach  9  shock  are  shown  in  Fig.  1(a).  The  CFD  solution  of  the  Navier-Stokes 
equations  predicts  a  shock  that  is  thinner  than  the  experimental  measurements  of  Alsmeyer  [8]  whereas  DSMC 
provides  almost  perfect  agreement.  The  hybrid  code  is  initialized  to  the  incorrect  CFD  result  and  “corrects”  it 
such  that  it  also  provides  excellent  agreement  with  the  measured  data.  Figure  1(b)  demonstrates  the  very  good 
agreement  obtained  with  the  hybrid  method  across  a  range  of  shock  Mach  numbers. 


FIGURE  1.  Argon  shock  waves:  (a)  normalized  density  profile  at  Mach  9;  (b)  reciprocal  shock  thickness. 

Figures  2  illustrate  the  performance  of  the  hybrid  code  for  two  dimensional,  Mach  12  flow  of  nitrogen  over  a 
cylinder.  In  Fig.  2(a),  significant  differences  for  the  temperature  contours  are  found  between  the  CFD  (lower) 
and  DSMC  (upper)  solutions.  Again,  the  hybrid  simulation  is  initialized  to  the  (presumably)  incorrect  CFD 
result,  and  then  modifies  it  to  provide  almost  perfect  agreement  with  the  full  DSMC  computation.  Indeed,  as 
illustrated  in  Fig.  2(b),  the  level  of  agreement  between  solutions  from  full  DSMC  and  the  hybrid  method  agree  at 
the  level  of  the  velocity  distribution  function.  Significantly,  in  the  best  case  we  have  found  thus  far  for  two- 


dimensional  configurations,  the  hybrid  code  produces  essentially  the  same  solution  as  full  DSMC  at  a  factor  of 
30  lower  overall  computational  cost  while  using  one  fifth  of  the  amount  of  computer  memory  [9].  Current  work 
is  focused  on  extension  of  the  hybrid  method  to  3D  and  further  research  will  be  required  to  extend  the  technique 
to  simulation  of  chemical  reactions. 


FIGURE  1.  Mach  12,  Kn=0.01  flow  of  nitrogen  about  a  cylinder:  (a)  contours  of  translational  temperature  obtained  with 
hybrid  (MPC),  DSMC,  and  CFD  methods;  (b)  velocity  distribution  functions  in  the  bow  shock. 

HYBRID  METHOD  FOR  NONEQUILIBRIUM  PLASMA  FLOWS 

Hall  thrusters  are  an  efficient  form  of  plasma  electric  propulsion  for  spacecraft.  Models  have  been  developed 
of  the  plasma  plume  of  such  thrusters  in  order  to  assess  spacecraft  integration  issues.  The  high  energy  ions 
created  by  the  thruster  can  sputter  spacecraft  surfaces  upon  impact  leading  to  possible  damage  and  subsequent  re¬ 
deposition.  In  this  section,  the  plume  of  a  Hall  thruster  is  modeled  using  a  hybrid  particle-continuum  approach. 

Model  Description 

Hall  thrusters  primarily  use  xenon  as  propellant.  The  plasma  at  the  exit  of  a  Hall  thruster  is  strongly 
nonequilibrium  with  different  species  possessing  different  temperatures  and  velocities.  Table  1  lists  typical 
thruster  exit  conditions  for  the  SPT-100  Hall  thruster  [10].  In  addition,  the  total  number  density  is  of  the  order  of 
1018  m'3  giving  a  Knudsen  number  greater  than  one  that  places  the  plasma  in  the  rarefied  flow  regime. 
Computational  analysis  of  Hall  thruster  plumes  is  regularly  performed  using  a  hybrid  particle-continuum 
formulation.  The  direct  simulation  Monte  Carlo  (DSMC)  method  [3]  models  the  collisions  of  the  heavy  particles 
(ions  and  atoms).  The  Particle  In  Cell  (PIC)  method  [11]  models  the  transport  of  the  ions  in  electric  fields. 
Overall,  a  hybrid  approach  is  employed  in  which  the  electrons  are  modeled  using  a  fluid  (continuum)  description. 

Plasma  Dynamics 

Hall  thruster  plume  models  employ  a  hybrid  approach  in  which  heavy  species  are  treated  using  particles  and 
the  electrons  are  considered  as  a  fluid.  Almost  all  previous  hybrid  models  reduce  the  electron  fluid  model  to  the 
Boltzmann  relation.  This  requires  that  the  electrons  be  collisionless,  currentless,  isothermal,  and  un-magnetized. 
All  of  these  assumptions  are  questionable  in  a  Hall  thruster  plume,  particularly  in  the  plume  near-field.  Despite 
the  simplicity  of  the  model,  these  hybrid  methods  have  been  quite  successful  in  simulating  the  far-field 
properties  of  a  number  of  different  Hall  thrusters.  As  mentioned  earlier,  the  ions  and  neutrals  are  treated  using  a 
combination  of  PIC  for  transporting  the  ions  in  electrostatic  fields,  and  DSMC  for  performing  collisions  and 
transporting  the  neutral  atoms. 

The  most  widely  used  fluid  electron  approach  for  plasma  plume  simulations  is  the  Boltzmann  model.  In  this 
approach,  quasi-neutrality  is  assumed,  which  allows  the  ion  density  to  represent  the  electron  density.  By  further 
assuming  that  the  electrons  are  isothermal,  collisionless,  and  un-magnetized,  and  that  their  pressure  obeys  the 
ideal  gas  law,  pe=nekTe ,  the  Boltzmann  relation  is  obtained  from  the  electron  momentum  equation: 
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where  ne  is  the  electron  number  density,  *  indicates  a  reference  state,  (j)  is  the  plasma  potential,  k  is  Boltzmann's 
constant,  Te  is  the  constant  electron  temperature,  and  e  is  the  electron  charge.  The  potential  is  then  differentiated 
spatially  to  obtain  the  electric  fields. 

In  the  present  work,  a  much  more  detailed  approach  is  employed  that  considers  all  three  steady  state 
conservation  equations.  The  electron  continuity  equation  is  written  as: 
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where  is  a  velocity  potential  and  the  right  hand  side  is  the  ionization  source  term  with  na  as  the  atomic  number 
density,  and  C,  the  ionization  rate  coefficient.  The  spatial  distribution  of  the  ion  particles,  treated  using  DSMC- 
PIC,  gives  the  electron  number  density,  ne ,  under  the  assumption  of  charge  neutrality.  This  allows  the  electron 
velocity  vector  to  be  determined  through  solution  of  Eq.  (3). 

The  electron  momentum  equation  is  written  as  a  modified  Ohm's  Law: 
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where  j  is  the  current  density  and  cr is  the  electrical  conductivity.  Equation  (4)  is  solved  for  using  the  charge 
continuity  condition: 

V.j  =0  (5) 

to  obtain  the  plasma  potential.  Finally,  the  electron  energy  equation  is  solved  to  obtain  the  electron  temperature: 
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where  Ke  is  the  electron  thermal  conductivity,  TH  is  the  heavy  species  temperature,  and  Ve  is  the  electron 
velocity. 


Collision  Dynamics 

There  are  two  basic  classes  of  collisions  that  are  important  in  Hall  thruster  plumes:  (1)  elastic  (momentum 
exchange);  and  (2)  charge  exchange.  Elastic  collisions  involve  only  exchange  of  momentum  between  the 
participating  particles.  For  the  systems  of  interest  here,  this  may  involve  atom-atom  or  atom-ion  collisions.  For 
atom— atom  collisions,  the  Variable  Hard  Sphere  (VHS)  [3]  collision  model  is  employed. 

Charge  exchange  concerns  the  transfer  of  one  or  more  electrons  between  an  atom  and  an  ion.  For  singly 
charged  ions,  the  following  cross  section  measured  by  Miller  et  al.  [12]  is  used: 

oCEX  (Xe,Xe+)  =  (-23.31og10(g)  + 142.2)0.8423  x  lCT20m2  (7) 

Also  reported  in  [9]  are  charge  exchange  cross  sections  for  the  interaction  where  a  doubly  charged  ion  captures 
two  electrons  from  an  atom.  These  cross  sections  are  less  than  a  factor  of  two  lower  than  the  values  for  the  singly 
charged  ions  at  corresponding  energies.  In  the  present  model,  it  is  assumed  that  there  is  no  transfer  of  momentum 
accompanying  the  transfer  of  the  electron(s).  This  assumption  is  based  on  the  premise  that  charge  exchange 
interactions  are  primarily  at  long  range.  It  is  further  assumed  that  atom-ion  momentum  exchange  cross  sections 
are  equal  to  the  charge  exchange  cross  section. 

Boundary  Conditions 

For  the  computation  of  Hall  thruster  plumes,  boundary  conditions  must  be  specified  at  several  locations:  (1) 
at  the  thruster  exit;  (2)  at  the  cathode  exit;  (3)  along  the  outer  edges  of  the  computational  domain;  and  (4)  along 
all  solid  surfaces  in  the  computational  domain. 

Several  macroscopic  properties  of  the  plasma  exiting  the  thruster  are  required  for  the  computations. 
Specifically,  the  plasma  potential,  the  electron  temperature,  and  for  each  of  the  particle  species  we  require  the 
number  density,  velocity,  and  temperature.  In  the  real  device,  these  properties  vary  radially  across  the  exit  plane. 
The  approach  to  determining  these  properties  involves  a  mixture  of  analysis  and  estimation.  The  basic 
performance  parameters  of  mass  flow  rate,  thrust,  and  total  ion  current  are  assumed  known.  The  neutrals  are 


assumed  to  exit  the  thrust  at  the  sonic  speed  corresponding  to  some  assumed  value  for  their  temperature.  Finally, 
divergence  angles  for  the  lower  and  upper  edges  of  the  exit  channel  must  be  assumed.  Combining  all  this 
information  then  allows  all  species  densities  and  the  ion  velocities  to  be  determined.  Determination  of  the 
properties  of  multiple  charge  states,  for  example  Xe2+  is  considered  in  the  present  study,  requires  knowledge  of 
the  current  fraction  of  that  state. 

Both  fluid  and  particle  boundary  conditions  are  required  at  the  outer  edges  of  the  computational  domain.  The 
usual  field  conditions  employed  simply  set  the  electric  fields  normal  to  the  boundary  edges  equal  to  zero. 
Similarly  the  gradients  in  electron  temperature  normal  to  the  surfaces  of  the  outer  boundaries  are  set  to  zero.  The 
particle  boundary  condition  is  to  simply  remove  from  the  computation  any  particle  crossing  the  domain  edge. 
The  results  presented  in  this  article  consider  the  thruster  to  be  operating  in  the  vacuum  of  space.  For  analysis  of 
thrusters  operated  in  a  vacuum  chamber,  the  finite  back  pressure  of  the  faciliy  is  included  by  simulating  a  fixed 
density  background  of  xenon  atoms  at  room  temperature.  These  atoms  can  collide  with  heavy  species  emitted  by 
the  thruster  and  so  can  make  an  important  contribution  to  the  charge  exchange  plasma.  In  addition,  the 
background  atoms  affect  the  neutral-electron  collision  frequency  that  appears  in  the  transport  coefficients  of  the 
electron  fluid  model. 

The  solid  wall  surfaces  of  the  Hall  thruster  are  also  included  in  the  computation.  Along  these  walls,  the 
plasma  potential  is  set  to  zero  and  zero  gradient  in  electron  temperature  is  employed.  Any  ions  colliding  with  the 
walls  are  neutralized.  Both  atoms  and  neutralized  ions  are  scattered  back  into  the  flow  field  from  the  surface  of 
the  thruster  wall  assuming  diffuse  reflection  at  a  temperature  of  300  K. 

Results 

For  illustration,  we  consider  the  plumes  generated  by  SPT-100  Hall  thrusters  employed  on  the  Russian 
Express  spacecraft.  A  variety  of  sensors  were  employed  on  board  the  two  spacecraft  to  characterize  the  effects  of 
firing  the  Hall  thrusters  on  the  spacecraft  operation  and  environment.  The  present  study  focuses  on 
measurements  of  ion  energy  distributions.  The  analysis  presented  briefly  here  is  described  in  more  detail  by 
Boyd  [10].  The  operating  conditions  of  the  SPT-100  Hall  thruster  considered  in  the  present  study  are  as  follows: 
flow  rate  =  5.3  mg/s,  discharge  current  =  4.5  A  (82%  is  attributed  to  ions  of  which  30%  by  current  fraction  is 
attributed  to  Xe2+),  discharge  voltage  =  300  V,  specific  impulse  =  1,600  sec.  The  computational  domain  extends 
more  than  10  m  axially  from  the  thruster  exit  and  10  m  radially  from  the  thruster  centerline  to  cover  all  of  the 
Express  probe  locations.  This  is  achieved  using  a  mesh  containing  190  by  175  non-uniform,  rectangular  cells. 
The  radial  mesh  spacing  at  the  thruster  exit  is  5  mm  giving  just  4  cells  across  the  exit  of  the  thruster.  In  a  typical 
computation  approximately  4  million  particles  are  employed  with  about  60%  representing  ions  (both  single  and 
double  charged).  The  neutral  atom  flow  is  first  allowed  to  reach  a  steady  state  by  using  a  large  time  step.  The 
ions  are  then  subsequently  introduced  with  a  time  step  of  about  10‘7s.  The  computations  reach  a  steady  state  for 
the  ions  after  about  15,000  iterations  and  solutions  are  then  averaged  over  a  further  10,000  iterations.  The  total 
computation  time  is  about  24  hours  on  a  personal  computer. 

In  Figs.  3(a)  and  3(b),  contours  are  shown  of  the  plasma  and  neutral  atom  number  densities,  respectively. 
These  show  that  the  two  populations  follow  quite  different  plume  expansion  dynamics.  The  bubble  of  charge 
exchange  plasma  formed  vertically  above  the  thruster  exit  plane  can  be  seen  clearly  in  Fig.  3(a). 

The  ion  energy  distribution  measured  on  Express  in  the  primary  beam  near  to  the  centerline  (at  7  deg.)  and  a 
distance  of  3.76  m  is  considered  in  Fig.  4(a),  where  the  Express  data  is  compared  with  the  results  of  the 
simulation.  Note,  in  terms  of  plotting  style,  that  exact  agreement  between  the  data  sets  would  mean  that  the  solid 
line  employed  for  the  model  results  would  go  through  the  center  of  the  horizontal  bar  of  each  column  of  the 
histogram  used  for  the  Express  data.  The  data  measured  in  space  provides  a  narrower  distribution  than  profiles 
measured  in  ground  facilities  and  this  is  perhaps  explained  by  collisional  broadening  present  in  the  vacuum  tank 
experiments.  In  Fig.  4(a),  there  is  clearly  very  good  agreement  between  the  simulation  and  the  Express  data. 

The  ion  energy  distribution  obtained  on  Express  at  the  large  angle  of  77  deg  and  a  distance  of  1.40  m  is  now 
considered.  This  location  is  of  particular  interest  since  it  is  characterized  primarily  by  charge  exchange  ions. 
Very  few  beam  ions  are  expected  to  exit  the  Hall  thruster  at  such  large  angles.  In  Fig.  4(b),  the  Express  data  are 
compared  with  the  results  from  the  simulation.  Figure  4(b)  includes  a  high  energy  structure  measured  on-board 
the  Express  spacecraft  that  extends  up  to  values  associated  with  primary  beam  ions  of  about  260  eV.  These  high 
energies  are  not  simulated  by  the  model,  although  the  peak  of  the  distribution  at  about  28  eV  is  very  well 
predicted. 
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FIGURE  3.  Plume  results  for  the  SPT-100  Hall  thruster:  (a)  Xe+  number  density;  (b)  Xe  number  density. 


FIGURE  4.  Ion  energy  distribution  functions  in  the  plume  of  a  Hall  thruster:  (a)  main  beam;  (b)  charge  exchange  plasma. 
Table  1:  Plasma  properties  at  the  exit  of  the  SPT-100  Hall  thruster. 


Property/Species 

Xe 

Xe+ 

Xe2+ 

E- 

Number  Density  (m'3) 

1.2xl018 

2.4xl017 

2.6xl016 

2.9xl017 

Temperature 

750  K 

1  eV 

1  eV 

6  eV 

Velocity 

280  m/s 

18  km/s 

25  km/s 

-6  km/s 

An  example  of  use  of  the  hybrid  method  to  analyze  a  Hall  thruster  experiment  performed  in  a  vacuum  facility 
is  shown  in  Fig.  5.  In  this  case,  a  high  power  thruster  operated  at  10  kW  is  in  the  12V  vacuum  chamber  is 
considered  [13].  Complex  geometric  structures  representing  baffles  and  pumping  surfaces  inside  the  chamber 
are  included  the  analysis.  Future  directions  in  the  development  of  this  technique  include  extension  to  3D, 
inclusion  of  magnetic  field  effects,  and  improved  consideration  of  the  thruster  cathode. 


121- 
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Faraday  Cup  Probes 


FIGURE  5.  Hybrid  analysis  of  a  10  kW  Hall  thruster  operated  in  a  vacuum  chamber  [13]:  (a)  mesh;  (b)  plasma  potential. 

SUMMARY 


Numerical  simulation  of  nonequilibrium  gas  and  plasma  dynamics  represents  a  multi-scale  problem  that  can 
be  addressed  with  physical  accuracy  and  numerical  efficiency  using  hybrid  particle-continuum  techniques.  For 
hypersonic  gas  dynamics,  a  hybrid  method  was  described  that  decomposes  a  flow  domain  into  regions  that  are 
simulated  using  continuum  (CFD)  or  particle  (DSMC)  methods.  This  technique  is  able  to  reproduce  results  from 
full  DSMC  computations  at  the  level  of  the  velocity  distribution  function  while  requiring  only  a  fraction  of  the 
cost.  For  plasma  jets,  a  hybrid  method  was  described  that  uses  a  particle  method  (DSMC-PIC)  for  the  heavy 
species,  and  a  fluid  approach  for  the  electrons.  This  technique  makes  it  possible  to  conduct  two-dimensional 
simulations  in  just  a  few  hours  that  offer  excellent  agreement  with  data  measured  during  space  operation. 
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Overview 


•  Nonequilibrium  gas  and  plasma  dynamics 

•  Modeling  of  hypersonic  aerothermodynamics: 

-  hybrid  fluid-particle  method 

-  illustrative  examples 

•  Modeling  of  spacecraft  plasma  propulsion  systems: 

-  effects  of  thruster  plumes 

-  hybrid  fluid-particle  method 

-  example:  Hall  thruster  plumes 


2 


Distribution  A:  Approved  for  public  release;  distribution  unlimited. 


What  is  Nonequilibrium? 


x(m) 


Continuum  flow  regime: 

-  small  Kn=A,/L 

-  high  density  and/or  large  L 

-  many  collisions 

-  thermodynamic  equilibrium 


x  (m) 

Rarefied  flow  regime: 

-  large  Kn=X/L 

-  low  density  and/or  small  L 

-  few  collisions 

-  nonequilibrium 
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Nonequilibrium  Systems 


Hypersonic  Vehicles 


Micro-scale  Gas  Flows 


Spacecraft  Propulsion 


Vacuum  Systems 
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Computation  of 
Nonequilibrium  Gas  Flow 


Flow 

continuum 

Regimes: 

◄  ◄ 

◄ 

Kn  0.01 

◄ 


slip  transitional  free-molecular 
►  ◄  ►  ► 

► 

0.1  1 
DSMC 


Model 

Accuracy: 


Boltzmann  Equation  Collisionless 


«  N-S 

Euler  Burnett  Boltzmann  Eqn 


Numerical  DSMC 

Cost: 

NS 

◄ 

Kn  0.01  0.1  1 
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Direct  Simulation 
Monte  Carlo  Method  (DSMC) 


Particle  method  for  nonequilibrium  gas  flows: 

-  developed  by  Bird  (1960’s) 

-  particles  move/collide  in  physical  space 

-  particles  possess  molecular  level 
properties,  e.g.  u’  (thermal  velocity) 

-  cell  size  Ax  ~  X,  time  step  At  ~  x 

-  collisions  handled  statistically  (not  MD) 

-  internal  energy  relaxation,  chemistry 

-  gas-surface  interactions 


► 


u\  v\  w’ 
X,  y,  Z 

^rot’  ^vib 
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Hypersonic  Vehicles: 

Ma  >  5 


Ballistic  Missiles 


Space  Planes 


Entry  Capsules 


Air-Breathing  Cruisers 
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Hybrid  Method 
for  Hypersonic  Flows 


•  Hypersonic  vehicles  encounter  a  variety  of  flow  regimes: 

-  continuum:  modeled  accurately  and  efficiently  with  CFD 

-  rarefied:  modeled  accurately  and  efficiently  with  DSMC 


Rarefied  DSMC  approach: 

based  on  kinetic  theory 
high  altitude,  low  density 
small  length  scales 

Continuum  CFD  approach: 

solve  NS  equations 
low  altitude 
long  length  scales 
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Hybrid  CFD-DSMC  Approach 


Scalabrin  and  Boyd 
AIAA-2006-3773 


Dietrich  and  Boyd 
J.  Comp.  Phys.,  126,  1996 


NS  solver  (LeMANS)  « 

-  2D/3D  unstructured  mesh 

-  Modified  Steger-Warming 
flux-vector  splitting,  parallel 

-  Point-implicit  time  integration 

-  Many  physical  models 


interface 


►  DSMC  (MONACO) 

-  2D/3D  unstructured  mesh 

-  Parallel,  domain  decomposition 

-  Many  physical  models 


► 


Pi?  Tr,  Tv 
U,  V,  w 


u\  v\  w’ 
X,  y,  Z 

^rot’  ^vib 
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Simulation  Codes 


•  CFD=LeMANS  code: 

-  developed  by  Scalabrin  &  Boyd  (2006) 

-  solves  the  Navier-Stokes  equations 

-  isothermal,  no-slip  wall  conditions  employed 

-  transport  coefficients  from  collision  integrals 


•  DSMC=MONACO  implementation: 

-  developed  by  Boyd  et  al.  since  1996 

-  general  2D/AXI/3D,  parallel,  unstructured  meshes 

-  isothermal  wall,  some  slip  always  present 

-  collision  model  consistent  with  CFD  transport 
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Continuum  Breakdown 


How  do  we  know  when  CFD  fails? 

-  continuum  breakdown  parameter 

-  gradient-length  Knudsen  number,  KnGL 

KnGL.Q=^Q\ 


KnGL  vai\\(Kn(i,_  n ,  KnGTV ,  KnGT_T ) 


" GL-p 


GL-V 


GL-T 


For  regions  in  the  flow  field  where  KnGL  <  0.05 
-  it  has  been  shown  that 


DSMCsol.n  -CFD 

dsmcS0L.n 


SOL'N 


<5  9, 
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Hybrid  Coupling  Method 


•  State-based  coupling 

-  Use  existing  CFD  and  DSMC  boundary  procedures 
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DSMC  Boundary  Conditions 


•  Generate  particles  in 
the  interior  and 
boundary  of  DSMC 
domain 

•  Use  the  Chapman- 
Enskog  distribution 
based  on  local  CFD 
information 


Velocity  profile  through 
a  normal  shock 
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CFD  Boundary  Conditions 


•  Average  DSMC  data 
communicated  to  CFD  has 
huge  fluctuations 


•  Sub-relaxation  technique 
(Sun  and  Boyd,  2005) 

A^d-ejA^  +  BAj 

•  History  before  time-step  i 
is  removed  by 


(1-0) 


J-i 


A A i  +  .  . 

3  3  l-(l-0);-'  3 

j=e+i 
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Hybrid  CFD-DSMC  Algorithm 


Initial  CFD  solution 
Decompose  DSMC  and  CFD  regions 
Generate  particles  in  DSMC  regions 


Iterate  DSMC  regions 
Use  sub-relax  avg.  to  track  macro  changes 

T 

Apply  KnGL  to  update  interface  locations 


T 

Update  CFD  BCs  and  converge 


o 


Have  interfaces  stopped  moving? 


Yes 


Lock  interfaces 
►  Sample  DSMC 
Converge  CFD 
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Example: 

Normal  Shock  Waves 


•  Argon  and  nitrogen  normal  shocks  investigated: 

-  relatively  simple  hypersonic  flow 

-  Alsmeyer  experimental  measurements 

-  well-known  case  for  testing  new  algorithms 

•  Simulations: 

-  modeled  in  2D  (400  x  5  cells) 

-  initialized  by  jump  conditions 

-  pure  DSMC 

-  pure  CFD  (Navier-Stokes  equations) 

-  hybrid  code  initialized  by  CFD  solution 
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Argon  Shock  Wave: 
Mach  9 


x/X 
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Physical  Accuracy 
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Y  [m] 


2D  Cylinder  Flow  of  N2 
Ma=12,  Kn,=0.01 


CFD  corrections  needed  in: 

-  bow  shock 

-  boundary  layer 
-wake 


iiS 

NG 

PD 

_aboratory 

|  Navier-Stoke 
■  DSMC 

5 

Hybrid  algorithm: 

-  initialized  by  CFD 

-  final  solution  agrees  with 
pure  DSMC 
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Local  Velocity  Distributions 


gNG 

2PD 

Laboratory 


Hybrid  method  is  physically  accurate  to  the  level  of  the 
particle  velocity  distribution 


Vx  [m/s] 


Vx  [m/s] 


VY[m/s] 


VY  [m  /  s] 


Mid-shock  along  stagnation  line 


1 35°  around  cylinder  surface 
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Future  Directions 
for  Hypersonic  Hybrid  Method 


•  Algorithm  development: 

-  generalized  mesh,  3D,  parallel  (Deschenes) 

-  DSMC:  implicit  and/or  other  acceleration  schemes 

-  other  approaches:  e.g.  LD-DSMC  (Burt) 

•  Physics  development: 

-  chemically  reacting  gas  mixture 

-  extension  to  plasma  flows  (MHD) 

•  Development  of  hybrid  methodology: 

-  refinement  of  breakdown  parameters 

-  evaluation  against  data  (measured,  computed) 
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Spacecraft 

Plasma  Propulsion  Systems 


Tasks:  orbit  transfer  and  maintenance 


Electrothermal  Propulsion 

Resistojet  thrusters 
Arciet  thrusters 

Electrostatic  Propulsion 

Ion  thrusters 
Hall  thrusters 


Electromagnetic  Propulsion 
Pulsed  plasma  thrusters 
Magnetoplasmadynamic  thrusters 
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Spacecraft  Integration: 
Effects  of  Thruster  Plumes 


•  Spacecraft-plume  interaction: 

-  ions  diverge  from  thruster 

-  energetic  impingement 

-  material  sputtering 

-  subsequent  re-deposition 

-  reduced  functionality 

•  Accurate  simulations: 

-  experiments  difficult 

-  Kn~1  at  thruster  exit 

-  charge  exchange 

-  electro-static  forces 
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Plume  Plasma  Dynamics: 
Particle  In  Cell  (PIC) 


•  Hybrid  method:  particles  for  ions/atoms  (nonequilibrium) 

fluid  for  electrons  due  to  small  mass 


Heavy  species  treated  using  PIC: 

-  weight-to-mesh  (Ruyten,  1993) 

-  charge  neutrality  =>  ne 

-  (|>  from  fluid  electron  models 

-  differentiate  (|)  for  electric  fields,  E 

-  weight-to-particles  to  accelerate 

-  neutral  transport  included 

-  combined  with  DSMC  for  collisions 


{u\  v\  w’ 
x,y,z 

m,  q 
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Plume  Electron  Dynamics: 
1.  Boltzmann  Relation 


Standard  approach  for  plasma  thruster  plumes  : 

-  derived  from  electron  momentum  equation 

-  currentless,  isothermal,  un-magnetized,  collisionless 

-  charge  neutrality  provides  potential  from  ion  density: 


kT , 

(j)  -  (f>  *  =  —  In 

e 

-  n*,  (|)*  are  reference  values 


n  x 


\n  */ 


Required  models/inputs: 

-  boundary  conditions  at  thruster  exit 
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Plume  Electron  Dynamics: 
2.  Detailed  Fluid  Model 


•  ADI  solution  of  steady-state  conservation  equations: 

-  electron  continuity  equation  (includinq  ionization) 

\  \p  =  nnaCi  where Wip  =  nu 


-  generalized  Ohm’s  law  (electron  momentum  equation) 


-  electron  energy  equation  (including  ionization) 

1  /  -  —  3  \ 

V  T  =  — Vlr^fO-VT^  H — i  - j.E  +  ~n(w.V)^71  +  pW.u  nnaCl£i j 

•  Required  models/inputs: 

-  transport/ionization  coefficients  for  propellant  system 

-  boundary  conditions  at  thruster  and  cathode  exits 
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Thruster  Exit  Conditions 


For  the  SPT-100 


Inner  Diameter 

60  mm 

Outer  Diameter 

100  mm 

Plasma  Density 

2.4xl017  m'3 

Neutral  Density 

2.0xl018m'3 

Ion  Velocity 

18.5  km/s 

Neutral  Velocity 

280  m/s 

Electron  Temperature 

6eV 

Ion  Temperature 

leV 

Neutral  Temperature 

750  K 
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R(m) 


Flow  Field: 

Heavy  Particle  Properties 


Ion  Number  Density  (nr3)  Atom  Number  Density  (nr3) 


28 


Distribution  A:  Approved  for  public  release;  distribution  unlimited. 


Ion  Energy  Distributions 


Primary  Beam 
(z=4  m,  0=7.5  deg.) 


CEX  Plasma 
(z=1 .4  m,  0=77.5  deg.) 
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Analysis  of 

The  12V  EP  Test  Facility 


gNG 

2PD 

Laboratory 


f 


1 


Left  -  30°  divergence,  Boltzmann 
Right  -  30°  divergence,  Detailed 
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Future  Directions 
for  Hybrid  Plume  Modeling 


•  Improved  numerical  algorithms: 

-  unstructured  3D  meshes 

•  Improved  physical  modeling: 

-  plasma  transport  coefficients 

-  thruster  exit  conditions  from  thruster  simulation 

-  effects  on  electrons  of  thruster  magnetic  field 

•  Further  assessments  of  the  models: 

-  other  Hall  thrusters  (e.g.  anode  layer  thrusters) 

-  other  EP  thrusters  (e.g.  MPD,  ion,  PPT) 
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Summary 


•  Hybrid  method  for  hypersonic  aerothermodynamics: 

-  combine  accuracy  of  DSMC,  efficiency  of  CFD 

-  decompose  domain  into  CFD  and  DSMC  regions 

-  hybrid  algorithm  successfully  “corrects”  CFD 

-  best  case  (so  far)  more  than  30  times  faster  than  DSMC 


•  Hybrid  method  for  spacecraft  thruster  plumes: 

-  impossible  to  model  electrons  as  particles  (mass) 

-  particles  for  ion/atoms,  fluid  for  electrons 

-  both  approaches  used  simultaneously  everywhere 

-  good  agreement  with  plume  measurements 
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